function cook_membrane
clc
%close all
% addpath('../');
% addpath('../../');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create mesh/read mesh
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% projDir = 'home/hbui/kratos_janosch';
% projName = 'cook-membrane.gid';
%fn = [projDir '/' projName '/' 'mesh4.msh'];

%set of quad mesh
% solve('/home/hbui/kratos_janosch/cook-membrane.gid/mesh4.msh')
% solve('/home/hbui/kratos_janosch/cook-membrane.gid/mesh9.msh')
% solve('/home/hbui/kratos_janosch/cook-membrane.gid/mesh16.msh')
% solve('/home/hbui/kratos_janosch/cook-membrane.gid/mesh25.msh')
% solve('/home/hbui/kratos_janosch/cook-membrane.gid/mesh36.msh')
% solve('/home/hbui/kratos_janosch/cook-membrane.gid/mesh49.msh')
% solve('/home/hbui/kratos_janosch/cook-membrane.gid/mesh64.msh')
% solve('/home/hbui/kratos_janosch/cook-membrane.gid/mesh81.msh')
% solve('/home/hbui/kratos_janosch/cook-membrane.gid/mesh100.msh')

%set of tri mesh
solve('/home/hbui/kratos_janosch/cook-membrane.gid/trimesh1.msh')

% rmpath('../');
% rmpath('../../');
end

function y = solve(fn)
mg = fem2.gid.GidMeshGenerator(fn);
mesh = mg.getMesh(2);
M = fem2.Model(mesh);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
meshUtil = fem2.MeshUtilities(mesh);
c1 = fem2.Constraint([0 0]); %fix x, fix y
leftedge = meshUtil.seekNodesOnSurface(1,0,0);
rightedge = meshUtil.seekNodesOnSurface(1,0,-48);
tip = meshUtil.seekNode([48 60]);
mesh.setConstraint(c1,leftedge);
F = 100;
for i=1:rightedge.size()-1
    ldx = math2.ConstantFunctionRnToR(0);
    ldy = math2.ConstantFunctionRnToR(F);
    n1 = rightedge.get(i-1);
    n2 = rightedge.get(i);
    e = mesh.addEdge(fem2.Edge([n1 n2]));
    el = fem2.EdgeLoad(e,[ldx ldy]);
    M.addLoad(el);
end
mesh.setNumberOfFields(2);
% plotMesh(mesh,{'normal','conditions','id'});
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% element declaration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nf = mesh.countFaces();
nu = 0.4999;
E = 250;
t = 1.0;
rho = 0.0;
ss = fem2.enu.State.PLANE_STRAIN;
cm = fem2.StVenantKirchhoffMM(E, nu, t, rho, ss);
for i=1:nf
    
    %option StructuralELement
%     e = fem2.StructuralElement(mesh.getFace(i-1),cm);
%     e.setIntegrationType(fem2.enu.IntegrationType.REDUCED_INTEGRATION);
%     e.setIntegrationType(fem2.enu.IntegrationType.SELECTIVE_REDUCED_INTEGRATION);

    %option BBarStructuralElement
    e = fem2.BBarStructuralElement(mesh.getFace(i-1),cm);
    
    M.addElement(e);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
M.preAnalysis();
M.Initialize();
n = mesh.countDofs();
kk = inf.jlinalg.Array2DMatrix(n,n);
M.assembleKt(kk);
r = zeros(1,n);
r = M.assembleR(r);
u = kk.getValues()\r;
M.setX(1,u);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% post process
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plotDeformedMesh(mesh,1);

y = tip.getUHatLocal()';
end
